GISS Assessment2

R Markdown

This is R Markdown file produced to evidence the spatial data analysis skills learned throughout CASA0005: Geographical Information System and Sience, specifically for part 2 of its assessment.

setwd("D:/Xin/MSc/Modules/CASA0005_GISS/Assessment/Assessment Part 2/")

library(geojsonio)

## 
## Attaching package: 'geojsonio'

## The following object is masked from 'package:base':
## 
##     pretty

library(tmap)
hunt <- geojson_read("https://www.dropbox.com/s/wa2ip35tcmt93g3/Team7.geojson?raw=1", what = "sp")
tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(hunt) + tm_lines(col="red", lwd=4)

## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3

##########################################
####Read in tube stations
library(sf)

stn <- st_read("https://www.doogal.co.uk/LondonStationsKML.ashx")

## Reading layer `London stations with zone information' from data source `https://www.doogal.co.uk/LondonStationsKML.ashx' using driver `KML'
## Simple feature collection with 641 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: -0.61143 ymin: 51.28216 xmax: 0.329952 ymax: 51.74702
## epsg (SRID):    4326
## proj4string:    +proj=longlat +datum=WGS84 +no_defs

ukgrid="+init=epsg:27700"
stn<-st_transform(stn,ukgrid)
tmap_mode("view")

## tmap mode set to interactive viewing

tm_shape(stn)+tm_markers(col="red",size=0.01,shape=10)

## Symbol shapes other than circles or icons are not supported in view mode.

##########################################
####geocode treasure locations

library(ggmap)

## Loading required package: ggplot2

## Google Maps API Terms of Service: http://developers.google.com/maps/terms.

## Please cite ggmap if you use it: see citation("ggmap") for details.

library(tidyverse)

## -- Attaching packages ---------------------------------------------------------------------------------- tidyverse 1.2.1 --

## v tibble  1.4.2     v purrr   0.2.5
## v tidyr   0.8.1     v dplyr   0.7.7
## v readr   1.1.1     v stringr 1.3.1
## v tibble  1.4.2     v forcats 0.3.0

## -- Conflicts ------------------------------------------------------------------------------------- tidyverse_conflicts() --
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(RCurl)

## Loading required package: bitops

## 
## Attaching package: 'RCurl'

## The following object is masked from 'package:tidyr':
## 
##     complete

library(RJSONIO)

## Warning: package 'RJSONIO' was built under R version 3.5.2

library(plyr)

## -------------------------------------------------------------------------

## You have loaded plyr after dplyr - this is likely to cause problems.
## If you need functions from both plyr and dplyr, please load plyr first, then dplyr:
## library(plyr); library(dplyr)

## -------------------------------------------------------------------------

## 
## Attaching package: 'plyr'

## The following objects are masked from 'package:dplyr':
## 
##     arrange, count, desc, failwith, id, mutate, rename, summarise,
##     summarize

## The following object is masked from 'package:purrr':
## 
##     compact

#Read in the list of treasure locations and their points scored from the dropbox given
huntaddresses <- read_csv("https://www.dropbox.com/s/v66l4cx7aia9jlo/Locations.csv?raw=1")

## Parsed with column specification:
## cols(
##   Location = col_character(),
##   Points = col_integer()
## )

#Set google map api for geocoding purposes
my_api="AIzaSyA5EUX3ZcX4u6b_OhogTM9FaWlxu4U5fJc"
register_google(key=my_api)
ggmap(get_googlemap())

## Source : https://maps.googleapis.com/maps/api/staticmap?center=29.763284,-95.363271&zoom=10&size=640x640&scale=2&maptype=terrain&key=AIzaSyA5EUX3ZcX4u6b_OhogTM9FaWlxu4U5fJc

#highlight this whole block and create this function to access the Google Places API
url <- function(address, return.call = "json", sensor = "false") {
  root <- "https://maps.google.com/maps/api/geocode/"
  u <- paste(root, return.call, "?address=", address, "&sensor=", sensor, "&key=", my_api, sep = "")
  return(URLencode(u))
}

#highlight this whole block and create this function to geocode some places just from a random list of treasure hunt locations
geoCode <- function(address,verbose=FALSE) {
  if(verbose) cat(address,"\n")
  u <- url(address)
  doc <- getURL(u)
  x <- fromJSON(doc,simplify = FALSE)
  if(x$status=="OK") {
    lat <- x$results[[1]]$geometry$location$lat
    lng <- x$results[[1]]$geometry$location$lng
    location_type <- x$results[[1]]$geometry$location_type
    formatted_address <- x$results[[1]]$formatted_address
    return(c(lat, lng, location_type, formatted_address))
    Sys.sleep(0.5)
  } else {
  return(c(NA,NA,NA, NA))
  }
}

#now use the geoCode() function (which calls the URL function) to geocode our list of places

#for loop to cycle through every treasure hunt location
i=1
for(i in 1:nrow(huntaddresses)){
  # Every nine records, pause 3 seconds so that the API doesn't kick us off...
  if(i %% 9 == 0) Sys.sleep(3)
  #now create a temporary list of useful elements
  tempdf <- as.list(geoCode(huntaddresses[i,1]))
  #and write these back into our dataframe
  huntaddresses[i,3] <- tempdf[1]
  huntaddresses[i,4] <- tempdf[2]
  huntaddresses[i,5] <- tempdf[4]
}
# rename the columns
names(huntaddresses) <- c("Location","Points","lat","lon","GoogleAddress")
head(huntaddresses)

## # A tibble: 6 x 5
##   Location            Points lat     lon    GoogleAddress                 
##   <chr>                <int> <chr>   <chr>  <chr>                         
## 1 Picadilly Circus, ~      2 51.509~ -0.13~ Piccadilly Circus, London W1B~
## 2 Nelson's Column, L~      2 51.507~ -0.12~ 5 Trafalgar Square, London WC~
## 3 Big Ben, London          2 51.500~ -0.12~ Westminster, London SW1A 0AA,~
## 4 100 Club, Oxford S~      5 51.516~ -0.13~ Century House, 100 Oxford St,~
## 5 Fabric, Charterhou~      5 51.519~ -0.10~ 77A Charterhouse St, Clerkenw~
## 6 The Gherkin, 30 St~      2 51.514~ -0.08~ 30 St Mary Axe, London EC3A 8~

#write a new .csv file to your working directory
write_csv(huntaddresses, "huntaddresses.csv")

##########################################
####Transform hunt data into simple feature objects, and measure the distance travelled on the hunt route

hunt_sf<-st_as_sf(hunt,ukgrid)
st_length(hunt_sf)

## 46610.68 [m]

#Transform hunt data into simple feature objects
hunt_sf<-st_transform(hunt_sf,ukgrid)

#Plot hunt route and stations on the same map
tm_shape(stn)+tm_markers(col="red",size=0.02,shape=16)+tm_shape(hunt_sf)+tm_lines(col="navy",lwd=4)

##########################################
####Create station buffers of 100m radii
stn_buff<-st_buffer(stn,dist=100)

#View station buffers and hunt route on the same map
tm_shape(stn_buff)+tm_polygons(col="red",size=0.02,shape=16,alpha=0.3)+tm_shape(hunt_sf)+tm_lines(col="navy",lwd=4)

#Count the instances of intersects between the station buffers created and the hunt route.
num_stn_intersects<-sum(lengths(st_intersects(stn_buff,hunt_sf)))
num_stn_intersects

## [1] 24

#Set latitude and longitude coordinate system
latlon<-"+init=epsg:4326"

#Retrieve coordinates data from huntaddresses data generated previously
huntaddresses$lat<-as.numeric(as.character(huntaddresses$lat))
huntaddresses$lon<-as.numeric(as.character(huntaddresses$lon))

#Create simple features from huntaddresses dataframe
huntadd_sf<-st_as_sf(x=huntaddresses,coords=c("lon","lat"),crs=latlon)

#And convert fro 4326 to 27700 crs
huntadd_sf<-st_transform(huntadd_sf,ukgrid)

#Check crs
st_crs(huntadd_sf)

## Coordinate Reference System:
##   EPSG: 27700 
##   proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 +units=m +no_defs"

#Create buffers for huntaddresses
huntadd_buff<-st_buffer(huntadd_sf,dist=300)

#View hunt location buffers and hunt route on one map
tm_shape(huntadd_buff)+tm_polygons(col="red",size=0.02,shape=16,alpha=0.3)+tm_shape(hunt_sf)+tm_lines(col="navy",lwd=4)

#Count the instances of intersects between the huntaddress buffers and the hunt route
num_huntadd_intersects<-sum(lengths(st_intersects(huntadd_buff,hunt_sf)))
num_huntadd_intersects

## [1] 16

##########################################
####Calculate max and min male life expectancy
#Read in the England shapefile as sp dataframe

library(geojsonio)
EnglandMap<-geojson_read("http://geoportal.statistics.gov.uk/datasets/8edafbe3276d4b56aec60991cbddda50_2.geojson", what = "sp")

#Crop the London shapefile only
LondonMap <- EnglandMap[grep("^E09",EnglandMap@data$lad15cd),]

#Quickly plot the London wards shapefile
qtm(LondonMap)

#Convert sp objects into simple features
LondonMap_sf<-st_as_sf(LondonMap)

#Convert the crs to 27700 and check
LondonMap_sf<-st_transform(LondonMap_sf,ukgrid)
st_crs(LondonMap_sf)

## Coordinate Reference System:
##   EPSG: 27700 
##   proj4string: "+proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.1502,0.247,0.8421,-20.4894 +units=m +no_defs"

#Read in male life expectancy data
library(readxl)
mle_raw<-data.frame(read_excel("leatbirthandatage65byukla201517.xls",sheet=2,col_names=TRUE,col_types=NULL,na="",skip=2))

library(tmaptools)

LondonMap_sf<-append_data(LondonMap_sf,mle_raw,key.shp="lad15cd",key.data="Area.Codes",ignore.duplicates=TRUE,ignore.na=TRUE)

## Data contains duplicated keys: NA

## mle_raw key variable "Area.Codes" contains NA's, which are ignored

## Over coverage: 411 out of 444 data records were not appended. Run over_coverage() to get the corresponding data row numbers and key values.

#Append indicator column for intersect
LondonMap_sf$IntersectIndicator<-0

#Assign 1 as indicator where wards intersect with hunt route
for(i in 1:nrow(LondonMap_sf)) {
  if(sum(lengths(st_intersects(LondonMap_sf[i,],hunt_sf)))==1) {
    LondonMap_sf[i,ncol(LondonMap_sf)]=1
  } else {
    LondonMap_sf[i,ncol(LondonMap_sf)]=0
  } 
}

#Slice the London wards where hunt route intersect
subset<-subset(LondonMap_sf,LondonMap_sf$IntersectIndicator==1)

#The intersecting London ward with the highest male life expectancy in 2015-2017 at births, and the years of expectancy
subset[which.min(subset$X2015.2017),]$lad15nm

## [1] Lambeth
## 380 Levels: Aberdeen City Aberdeenshire Adur Allerdale ... York

max(subset$X2015.2017,na.rm=T)

## [1] 82.66597

#The intersecting London ward with the lowest male life expectancy in 2015-2017 at births, and the years of expectancy
subset[which.max(subset$X2015.2017),]$lad15nm

## [1] Westminster
## 380 Levels: Aberdeen City Aberdeenshire Adur Allerdale ... York

min(subset$X2015.2017,na.rm=T)

## [1] 78.69568

#Calculate the average years of male life expectancies in 2015-2017 at births
sum(subset$X2015.2017,na.rm = T)/nrow(subset)

## [1] 71.08626

##########################################
####Test spatial randomness of hunt locations
#create London boundaries with Northing/Easting, and set a window as the London boundary

library(spatstat)

## Warning: package 'spatstat' was built under R version 3.5.2

## Loading required package: spatstat.data

## Warning: package 'spatstat.data' was built under R version 3.5.2

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## Loading required package: rpart

## 
## spatstat 1.57-1       (nickname: 'Cartoon Physics') 
## For an introduction to spatstat, type 'beginner'

library(rgdal)

## Warning: package 'rgdal' was built under R version 3.5.2

## Loading required package: sp

## rgdal: version: 1.3-6, (SVN revision 773)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 2.2.3, released 2017/11/20
##  Path to GDAL shared files: C:/Users/Xin/Documents/R/win-library/3.5/sf/gdal
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 4.9.3, 15 August 2016, [PJ_VERSION: 493]
##  Path to PROJ.4 shared files: C:/Users/Xin/Documents/R/win-library/3.5/sf/proj
##  Linking to sp version: 1.3-1

LondonMapBNG<-spTransform(LondonMap,ukgrid)
summary(LondonMapBNG)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##        min      max
## x 503576.3 561958.7
## y 155850.7 200934.0
## Is projected: TRUE 
## proj4string :
## [+init=epsg:27700 +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717
## +x_0=400000 +y_0=-100000 +datum=OSGB36 +units=m +no_defs
## +ellps=airy
## +towgs84=446.448,-125.157,542.060,0.1502,0.2470,0.8421,-20.4894]
## Data attributes:
##       lad15cd                   lad15nm            lad15nmw     objectid  
##  E09000001: 1   Barking and Dagenham: 1                :33   Min.   :294  
##  E09000002: 1   Barnet              : 1   Abertawe     : 0   1st Qu.:302  
##  E09000003: 1   Bexley              : 1   Blaenau Gwent: 0   Median :310  
##  E09000004: 1   Brent               : 1   Bro Morgannwg: 0   Mean   :310  
##  E09000005: 1   Bromley             : 1   Caerdydd     : 0   3rd Qu.:318  
##  E09000006: 1   Camden              : 1   Caerffili    : 0   Max.   :326  
##  (Other)  :27   (Other)             :27   (Other)      : 0                
##  st_lengthshape   st_areashape      
##  Min.   : 8929   Min.   :  2897649  
##  1st Qu.:28384   1st Qu.: 26797942  
##  Median :37664   Median : 37628571  
##  Mean   :39255   Mean   : 47682317  
##  3rd Qu.:46679   3rd Qu.: 56413925  
##  Max.   :74641   Max.   :150125298  
## 

window<-as.owin(LondonMapBNG)

#Create a point pattern object for the hunt locations
huntadd.ppp<-ppp(x=st_coordinates(huntadd_sf)[,1],y=st_coordinates(huntadd_sf)[,2],window=window)

## Warning in ppp(x = st_coordinates(huntadd_sf)[, 1], y =
## st_coordinates(huntadd_sf)[, : 1 out of 50 points had NA or NaN coordinate
## values, and was discarded

## Warning: 2 points were rejected as lying outside the specified window

huntadd.ppp<-as.ppp(huntadd.ppp)
plot(huntadd.ppp,pch=16,cex=1,main="Hunt Locations in London")

#Plot the kernel density estimation map
plot(density(huntadd.ppp, sigma = 1000))

#Run a simple Ripley's K test on the hunt location data in London boundar, the maximum radius plotted is set at 5km.
K<-Kest(huntadd.ppp,correction="border",rmax=5000)
plot(K)